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AN ENGINEERING OPTIMIZATION METHOD WITH APPLICATION TO 


STOL-AIRCRAFT APPROACH AND LANDING TRAJECTORIES 

Heinrich G. Jacob* 

Ames Research Center 


SUMMARY 


An optimization method has been developed that computes the optimal open loop inputs for a 
dynamical system by observing only its output. The method reduces to static optimization by 
expressing the inputs as series of functions with parameters to be optimized. Since the method is 
not concerned with the details of the dynamical system to be optimized, it works for both linear 
and nonlinear systems. This report explains the method and applies it to optimizing longitudinal 
landing paths for a STOL aircraft with an augmented wing. Noise, fuel, time, and path deviation 
minimizations are considered with and without angle of attack, acceleration excursion, flight path, 
endpoint, and other constraints. The method is easy to learn and easy to use. 


INTRODUCTION 


This report presents the results of a study of approximately optimum longitudinal landing 
trajectories for a short takeoff and landing (STOL) aircraft. Since the study was motivated by the 
desire to understand the character of the trajectories under various criteria, the approximate 
optimization procedure discussed in this report appeared appropriate. 

There are numerous methods for determining optimum open-loop trajectories. Generally, 
approaches that reduce optimization to a programing problem by parameterizing variables 
parameterize the states as well as the controls. The method employed in this study, however, 
parameterizes only the input controls. While the results obtained in this way are also approximate, 
they are obtained with a minimum of complexity. The approximate method employed in this 
study, then, is more of an engineering approach. It is neither elegant nor rigorous, but it is easy to 
use. While it may require more trials on a digital computer than more sophisticated methods, the 
complete setup, from initiation of the study to its completion, should be considerably shorter. 

The procedure is described in two sections. The first describes the operation of the 
optimization algorithm itself, how to choose directions and step sizes for search, how to optimize 
along each direction, and how constraints are to be handled. The second section shows what must 
be done to use the method in the landing problem. It discusses the mathematical model of the 
aircraft, the performance criteria, the three control variables, and, most importantly, the selection 
of the set of finite approximations to the optimum control function. 


♦National Research Council Postdoctoral Research Associate in residence at Ames Research Center. 



The final section of the report presents the results of the investigation. The first results 
concern straight line flight paths, minimum time, minimum fuel, and minimum noise trajectories. 
Finally, more practical results are obtained after the introduction of supplementary constraints 
associated with passenger comfort and pilot preference. 


THE OPTIMIZATION METHOD 


Before the particular problem of optimizing the flight paths is treated, the method of 
optimization will be described. First, the general idea of the approach will be discussed. Then the 
particular stages of the process will be considered in a little more detail. The optimization algorithm 
used in this study will be explained briefly with a more complete description given in an appendix. 


Principle of Operation 


The usual problem in optimizing the performance of a dynamic system is to find the control 
function time histories that drive the system so as to minimize some quality criterion. Here the 
problem is transformed to a static one of parameter optimization by expressing each control 
function as a finite series of known time functions. The optimal parameters of the series are then 
determined by an iterative optimization procedure. The type of series and the number of its terms 
are chosen so that the expected optimal response can be approximated with sufficient precision. 


Figure 1 shows the control functions, represented by a set of coefficients, Cjj(j) driving the 
dynamic system. The subscript j(i) indicates that the numbers of coefficients of the different 


Constraints 



Figure 1.— The optimization procedure. 


control variables need not be the same. The 
performance index, PI, evaluates the output of 
the system in terms of the chosen quality 
criterion. The evaluation is stored and compared 
with values previously determined. This 
comparison triggers a new selection of 
coefficients in a particular way, and the process 
is repeated until the performance index ceases to 
change. Constraints may be considered either by 
adding penalty functions to the quality criterion 
or by introducing boundaries directly into the 
search space over which the nonlinear 
programing algorithm operates. 


The type of series expansion chosen to represent the control function is conditioned by any a 
priori information available on the shape of the optimum functions. If no prior information is 
available, the simplest series representation to use is the polynomial approximation. 


u(t) = Cj + c 2 t + • • • + c n t n 1 


( 1 ) 


2 






However, because the accuracy of this representation is not uniform in the argument space (ref. 1 ), 
such other expansions as trigonometric functions or Legendre polynomials, whose orthogonality 
properties speed convergence, can also be used. A sine expansion over an interval [0, Tp] takes the 
form 

u(t) = Cj sin (t7r/Tp) + • • • c n sin (nttf/Tp) 

To save the computer time required to determine the sine for all n values, equation (2) is written in 
a trigonometric power series in R = ttf/Tp 

u(t) = c, [sin (R)] + c 2 [2 sin (R) cos (R)J 
+ c 3 [3 sin (R) - 4 sin 3 (R)] 

+ • • • 

+ c n |(n-l) sin (R) cos n_2 (R) - {(n-l.)/3] sin 3 (R) cos n_4 (R) 

+ [(n — 1 )/5] sin 5 (R) cos n_6 (R) -••• + •••} (3) 

Any information available on the shape of the optimal solutions can be used in selecting the 
approximating functions. An example of useful information is knowledge that the solution switches 
from one fixed level to another. Then the best choise of coefficients to be found would be the times 
of occurrence of the levels of control. The choice of switching times as parameters in this case saves 
computer time and can yield the exact optimal control function (ref. 2). 

The performance index is the number computed by evaluation of system output in terms of 
the quality criterion. Since the number is calculated at each trial of each stage in optimization, the 
computation should be made rapidly. 

Constraints (either equalities or inequalities) may be imposed on the input, output, or internal 
state variables of the system. If the constraints are inequalities, transition from a permissible to a 
forbidden region may be gradual or abrupt. 

The method of optimization used for this study can handle both types of constraints. For 
example, the parameters of the control functions are specified at each trial and applied to the 
system. Each time an inequality constraint violation is detected, the optimization algorithm is 
signaled to provide a new set of parameters until a set is obtained that violates no boundaries. 

Equality constraints, such as terminal constraints on the output, are handled by adding to the 
quality criterion penalty functions involving the conditions. Since the optimization algorithm 
involves quadratic functions in an interpolation scheme, the use of quadratic penalty functions is 
particularly convenient. Penalty functions, of course, may also be used for inequality constraints, 
especially those where the boundaries are approached gradually. 

Constraints imposed on the input functions can sometimes be handled directly by an 
appropriate choice of type of function for an input. For instance, control functions with switching 
time as a parameter can easily account for constraints on the number of switches allowed or their 
minimum time separation. 


3 



The Optimization Algorithm 


The value of the quality criterion, which is the index of performance of the system forgiven 
conditions, is a real valued function of several variables, namely, the parameters of the control 
function. The object of the optimization algorithm is to start with a given choice of parameters and 
from there to find the parameters that give an extreme (maximum or minimum) value to the 
criterion (fig. 2). Several algorithms are available to do this: the one selected should (1) converge 



C| 


Figure 2.— Performance index as a function of two 
coefficients. 


quickly to the nearest relative extreme point 
insensitive to curved valleys and sharp ridges in 
the parameter-criterion space; (2) permit the 
imposition of various types of constraints; and 
(3) be as simple and small as possible to save 
time for the user and the computer. The 
algorithm used in this study has these 
properties: It is simple and short; it appears to 
be less sensitive to ridges and valleys than other 
static nongradient methods; and it can handle 
inequality constraints on the control parameters 
or their functions. (The Fortran IV program for 



the algorithm is given in the appendix.) 

Figures 3 and 4 illustrate the operation of 
the algorithm in solving the problem of how to 
change the present or current values of the 
parameters to improve the performance index. 
Figure 3 shows a current choice of coefficients 
C c as a point in a plane of two parameters Cj 
and C 2 . This point has been reached from a 
previous point by a change in the parameter 


Figure 3.— Determination of the main search direction for 
a two coefficient problem. 


values along the line labeled “main direction of 
curre nt stag e.” Before another change is made to 
point C c+2 along the main direction of the next 
stage, the next main direction must be found by 



a number of changes along secondary directions 
as part of the current stage of computation. The 
number of secondary directions is one less than 
the number of coefficients to be determined. 

The point C c +j from which the first 
secondary direction will be drawn is found by 
changing the coefficients a predetermined 
amount DC first parallel to the main direction of 
the current stage, then antiparallel to this 


Figure 4.— Extrapolation along a line. 


direction. The performance index associated 


with each of the three points C c — DC, C c 


C c + DC is calculated as illustrated in figure 4, and a parabola is fitted to these three points. At 
point C c +! in figure 4 the interpolating (or extrapolating) parabola reaches its extreme value. This is 
the starting point for the first secondary direction, which is orthogonal to the main direction C c , 
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C c+ ; as a simple way of guaranteeing that all independent directions in parameter space are tested. 
This process of finding a starting point and a secondary direction is repeated until the number of 
main and secondary directions equals the number of (independent) coefficients o f the control 
functions. The final step of this stage of the calculation is to find the point in figure 3 C c+2 which is 
the point of the maximum of the parabola fitted through the performance ind ices o f the three 
points along the last secondary direction. The direction defined by points C c and C c + 2 is the main 
direction for the next stage of the calculation. 

The step size DC is variable during the computations. If any new interpolated point along a 
line is closer to the old point than one quarter of the current step size along this line, then the 
current step size will be divided by 4. On the other hand, if the new point is calculated to be farther 
away from the old point than 20 times the current step size, the step size will be multiplied by 2. 

Constraints limiting the evolution of the input, state, and output variables are taken into 
consideration as follows: Before and during any computation of the performance index, the 
variables are checked continually to make sure that they lie in the permitted area. As soon as one or 
several variables violate a boundary, the computer-control flow is given back to the optimization 
algorithm, which delivers another set of input function coefficients. This process is repeated if 
necessary until input functions are delivered that do not violate these boundaries. 

By means of the search along the continuously corrected optimal direction, the extrapolation 
or interpolation procedure, and an automatic adjustment of the variable step size DC along each 
direction, the algorithm converges quickly to the optimum and is able to follow accurately sharp 
ridges or steep-curved valleys in the coefficient space. (Further details are provided in the 
appendix.) 

From this general description, the dynamic optimization scheme can be seen to have the 
following characteristics. It is a direct method because the cost function is optimized by a direct 
evaluation of the corresponding quality criterion to index the performance. The optimization is 
open loop in the sense that each new set of initial states of the dynamical system requires a new 
determination of the optimal control function coefficients. The optimal controls form only an 
approximate solution to the problem since the 
number of coefficients is limited and the type of 
series expansion employed may not be the most 
succinct expression of control possible. The 
optimization procedure, furthermore, reveals only 
the local optimum nearest the initial set of 
coefficients. Broader claims for the optimum will 
necessitate the use of other control function 
expansions with new coefficients. 

ELEMENTS IN FLIGHT PATH OPTIMIZATION 

The Aircraft 

The augmentor wing STOL vehicle (fig. 5) 
considered in this study is a modified version of a Figure 5.- C-8A augmentor wing STOL aircraft. 
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Figure 6.— Augmentor flap arrangement. 


C-8 Buffalo aircraft by deHavilland of Canada 
(ref. 3). The wings and engine were modified to 
enhance the aircraft’s short takeoff and landing 
capability. The normal lift of the wings was 
augmented by the addition of a cold air duct 
and redesign of the flaps. The flaps are in effect 
a diffuser, and cold air ducted along the wing is 
exhausted through long nozzles (fig. 6). As the 
flow of air passes through the diffuser formed 
by the flap, it induces additional flow through 
the intake portion of the flap. The cold air flow 
effects added direct thrust from the throat of 
the flap and added lift due to increased 
aerodynamic circulation (ref. 4). 


Two Rolls Royce Spey Mk 801 SF jet engines replace the previous turboprop engines. The hot 
exhaust gas of the engines can be deflected through a total angle of 98° to provide additional lift 
and braking for the vehicle. The hot thrust vectoring and augmented'lift give this vehicle an unusual 
control capability. 


Mathematical Models for Computation 

The following simple two-dimensional equations of motion were used to model the aircraft. 

T 

it = Xa — g sin 6 + — cos a - 0w (4a) 

A m 

T 

w = Z a + g cos 6 sin a + du 

m 

x = u cos d + w sin d 
z = -u sin d + w cos d 

where, with the coordinates defined in figure 7, 
it acceleration along the body x-axis 

w acceleration along the body z-axis 

x horizontal geographical speed 

z vertical geographical speed 

g gravity acceleration 32.1725 ft/sec 2 

(9.805 m/sec 2 ) Figure 7 - Coordinate system. 



(4b) 

(4c) 

(4d) 
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6 pitch angle 

o thrust angle 

T hot thrust, lb (N) 

m W/g aircraft mass with W = 40,000 lb (17,438 kg) 

The forces and are given in terms of the lift and drag coefficients Cl and Cp> and the 
parameters and variables indicated below. 

X A = “CdC^ Cj, 5p) q — - C - ^ - a + C L (a, Cj, 5p) q S ^ “ (5a) 

Z A = Cj, 5 f ) q ■ S - - ^ - C L (a, Cj, 5p) q — ^ — (5b) 

where 

q = | P V2 

and 

p atmospheric density 0.00238 lb sec 2 /ft 4 (1.231 kg/m 3 ) 

V total speed, ft/sec (m/sec) 

S reference wing area 865 ft 2 (80.36 m 2 ) 
a angle of attack 

7 flight-path angle 

5p flap angle 

Cp Cl(«, Cj, 5p), coefficient of lift 

Cp Cj)(a, Cj, 5p), coefficient of drag 

Cj T^CO/qS, cold thrust coefficient 
T^(T) cold thrust, lb, a function of the hot thrust 

The control variables in equations (4) are: pitch angle 9; hot exhaust or thrust angle a; and throttle 
angle 0 e . The study assumed that the pilot or autopilot could generate 6 directly, with no 
intervening dynamics. The throttle angle 6> e is related to the hot thrust by 


7 



Thrust, lb (N) 


T = 355.1° 0 e (= 1592 0 e , N) 


( 6 ) 


with 0 e given in degrees. The flap setting 5p is taken to be constant at 75° for this study. 


Normal 


Normal 


With the above equations all state and output variables of the aircraft, necessary for the 
determination of the performance index, may be computed, with the exception of the parameters 
T(^(T), Cp (a, Cj, 5p) and C£)(a, Cj, 5p) whose values can be obtained from functions fitted to test 

data. Figure 8 (taken from Boeing Document 
D6-26065 TN, 1971) shows the hot and cold 
thrust as a function of the throttle angle. The 
aerodynamic lift and drag parameters Cp and 
Cj) are shown as a function of angle of attack a 
and cold thrust coefficient Cj in figure 9; the 
measured data (denoted by A) are also from 
Boeing Document D6-26065 TN. The cold thrust 
Tp is given empirically by the following expressions 



T c = T(1 - 0.0143 0 e ) 


= T(0.2557 + 8.57/0 e ) 


0 e < 20° (7a) 


0 e > 20° 


(7b) 


angle. 




(a) Aerodynamic lift. 


(b) Aerodynamic drag. 


Figure 9.— Aerodynamic lift and drag parameters for flap angle 6p = 75° and aileron angle 6^ = 0°. 
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The parameter identification method described in reference 5 was used to give the functional fit 
shown in figure 9 for Cl and Cp. The functions used for the fit were hyperbolas and polynomials. 
When coupled with some least squares criterion program, the optimization algorithm described in 
the appendix could also have been used to fit the data. 

After the data shown in figure 9 were obtained, the aircraft design was modified so that the 
ailerons droop automatically as the flaps are deflected. With the flaps set at 75°, the ailerons droop 
30°. This droop has a negligible effect on the drag coefficient, but changes the lift coefficient in the 
following way: 


Cl = Cl* + 0.1 15 + 0.28 Cj (8) 

where Cl anc * Cl* are the coefficients with and without droop considered, respectively, for a flap 
setting of 75°. 

A number of control and state constraints were imposed on the problem of determining 
optimum landing flight paths. Some are physical or hardware constraints, referred to as design 
constraints, that were taken into account in every calculation: throttle angle, thrust angle, angle of 
attack, and total flight speed. These constraints are listed in table 1 . 

Other constraints, referred to as safety and comfort constraints, are concerned with pilot’s 
preference or ride quality (table 2). For some computations of the flight trajectories the controls 
were changed only at certain times and remained constant otherwise, leaving the pilots free to check 
other flight parameters. These restrictions on the handling capabilities of the input commands were 
accounted for by the choice of an appropriate class of possible input functions. 


TABLE 1.- DESIGN 

CONSTRAINTS 

TABLE 2.- SAFETY AND 

COMFORT 

Throttle angle, d e 

0° < d c < 35° 

CONSTRAINTS 


Thrust angle, a 

18° < a < 116° 

Normal acceleration change 

|Aa n | < 0.1 5g 

a 

|d| < 60° /sec 

Pitch rate 

|0| < 10° /sec 

Angle of attack, a 

-10° 30° 

Throttle position 

16° < 0 e < 35° 

Forward speed, V 

V < 95 knots 

Angle of attack 

a < 8° (except 
during flare) 


Two different structures were chosen to represent the control variables. The first structure, 
allowing general unconstrained flight paths, was a modified finite sine series of the form: 

6 = c(l,l) + c(l,2) 4 I +c(l,3)S R + c(l,4)2S R C R + c( 1,5)(3S R - 4S r 3 ) 


+ c(1,6)(8S r C r 3 — 4 S r C r ) + c(1,7)(16S r C r 4 - 12S R C R 2 + S R ) 
+ c(l,8)(6S R C R 5 — 20 S r 3 C r 3 + 6 S r s C r ) 


9 



a = c(2,l) + c(2,2) 


x 


+ • • 


(9b) 


-Xi 

6q = c(3,l) + c(3,2) — — + • • • (9c) 

-xi 

where c(ij), i = 1 , 3; j = 1,8, coefficients whose optimal values have to be determined 
x geographical horizontal distance in feet (meters) to the touchdown point 

Xj initial geographical horizontal distance in feet (meters) to the touchdown point 
Sr sin(-7rx/xi) 

Cr COS (-7rx/xj) 

Formulas (9a), (9b), and (9c) are derived from equation (3) by adding a constant term and a linear 
term (ramp), and choosing the geographical horizontal distance x to the touchdown point as the 
independent variable. 



Figure 10.- Structure of input control functions 
comprising hyperbolic tangent. 


As will be seen, the optimal control 
histories resulting from the sine series for input 
functions present continuously varying and 
sometimes oscillating time histories. These are 
difficult for pilots. To obtain control functions 
acceptable to pilots the following structure of 
input functions was selected, which yields 
control time histories that vary only periodically 
(fig- 10). 


c( 1 , 1 ) — c( 1 ,4) t a nh 2 1 x — c( 1 ,3 ) j + c(l,l) + c(l,4) ^ c(l,4) - cQ,?) tjnh 2 [x - c(l,5)j 


Ic(l, 2)1 


lc(l, 6)1 


c(l,4) — c(l,7) c(l,7) — c( 1 , 1 0) , . 2[x — c(l,9)] c(l,7)-c(l,10) 

+ tanh _ 

2 2 lc(l, 8)| 2 


(10a) 


c(2,l) — c(2,4) 


(10b) 


c(3,l) — c(3,4) 
2 


(10c) 
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The 10 coefficients of the input control function shown in figure 10 permit three successive 
flight-path changes. As used in this study, each of the three control functions can have the same 
amount of freedom, as shown, with the one limitation that all three functions must switch at the 
same time. If fewer flight-path changes are foreseen, only seven or even four coefficients are 
necessary for each input control function. 


Definition of the Quality Criteria 

The following quality criteria were considered for the optimization. 

Time minimum— The corresponding performance index for any flight path considered requires 
simply measuring the time needed to move from the fixed initial point to the desired endpoint. 




dt = Tp 


(ID 


Fuel minimum— Vox this criterion it was assumed that the hot thrust is approximately 
proportional to the fuel rate. 

• tf 


PI = Cp / T(t)dt = Fp 


( 12 ) 


where 

T hot thrust in lb (N) 

tptf initial and final flight time; tf free 

Cp 0.22X 10~ 3 lb/lb sec (conversion factor) 

Noise minimum - From the numerous and complex noise definitions available, an approximate 
formula was developed to estimate the amount of noise inconvenience encountered by the people 
living in the community surrounding a STOL airport and by the people in the neighborhood of the 
runway. The performance index for the minimum noise criterion involves three major factors: the 
distance function, the expression used for maximum perceived noise, and the modification of the 
maximum perceived noise expression by which it is converted to the maximum effective perceived 
noise expression. 

In calculating noise, two expressions are used for the distance: the altitude of the airplane 
from the initial point until the airplane is 1550 ft (472 m) horizontally from the touchdown point; 
and the slant distance from the airplane to the ground 200 ft (61 m) lateral to the airplane’s ground 
track from the 1550-ft (472 m) point until touchdown. Two expressions are used because 
STOLports (refs. 6, 7) are considered to have a cleared area beginning at a point about 1550 ft 
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Figure 1 1 Noise level estimation. 


(472 m) from touchdown; touchdown is 
assumed to occur 550 ft (168 m) from the 
beginning of the runway pavement, and the 
cleared area begins 1000 ft (305 m) from the 
pavement. The noise is evaluated along a line 
200 ft (61 m) away from the centerline of the 
cleared area. In figure 1 1 a trajectory in the x-z 
plane is shown on the left and the line of noise 
evaluation in the x-y plane on the right. The 
population density is assumed uniform along the 
line of evaluation and the same for both distance 
expressions. 


The maximum perceived noise level (MPNdB) may be expressed approximately by (ref. 8): 


MPNdB = 123 + 25 1og lo (500/D N )+521og 10 (T/12,200) (13) 

where 123 is the perceived noise level in dB encountered at a distance of Djq = 500 ft (152 m) from 
the STOL aircraft flying with a hot thrust of T= 12,200 lb (54,680 N). 


The value of MPNdB varies along the flight path with the hot thrust T and the altitude of the 
aircraft. The distance Djq is measured from the noise source (aircraft) to the closest point along the 
flight path where people might be. From the previous discussion, this closest distance is taken to be 


D N = -z 

D N =</z 2 +200 2 


for Xj < x < -1 550 (-472) 
for —1550 < x < 0 


(v/z 2 + (61 ) 2 ) (—472 < x < 0) 


where x is the (negative) distance to the desired touchdown point and z is the (negative) altitude, 
both in ft (m) (fig. 1 1). The total integral amount of MPNdB along the complete flight path would 
be 

* f 

/ MPNdB(x)dx = J MPNdB(t)V x (t)dt(dB ft) (dB m) (14) 

xj h 

where V x = horizontal speed in ft/sec (m/sec). A time term is added to equation (14) to obtain the 
duration of the effective perceived noise. The effect of the noise duration on the sensation of 
people can be approximated by adding 3 to 8 dB (depending on the type of the noise) to the basic 
PNdB value for each doubling of the noise duration (refs. 9, 10) which yields 


EPNdB ~ PNdB + k log, 0 (t/t ref ) 


(15) 


with 10 < k < 27. In consequence, from expressions (13), (14), and (15) and knowing that the time 
is inversely proportional to the horizontal speed of the aircraft, one gets for the “integral maximum 
effective perceived noise” (/MEPNdB-dx) over the total flight path: 
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Np = 


-f 


x f 


MEPNdB(x)dx 


x i 


■n 


MPNdB(t) + k log 


■0 


-V x (t)dt (16a) 


Choosing k = 20 and taking as a reference speed V re f = 60 knots « 101 .2 ft/sec (30.85 m/sec) gives 
the index of performance corresponding to a minimum noise quality criterion 


-/ 


PI = I MEPNdB(t)V x (t)dt 


■a 

h 


MPNdB(t) + 20 log j o 


101.2 

V x (t)_ 


V x (t)dt = Np 


(16b) 


This minimized noise performance index Np has a dimension EPNdB feet (EPNdB meters) that 
is not very revealing. Hence, the value Np was divided by the horizontal flight distance to get an 
average effective noise Np encountered on ground from the initial (xj) to the final (xf) flight point: 

Np 

Np = (16c) 

x f- x i 

Deviation minimum— This quality criterion is based on the measurement of the integral 
absolute error between an actual flight path and a predefined, desired one. The desired path may 
consist of two straight line parts, for example, the first part of the trajectory might have a flight 
path angle of 7 = 0 ° (horizontal flight) and the second part (descent to the aim point), with a flight 
path angle of 7 = —7.5°. The index of performance corresponding to this criterion is calculated by 

PI = f |z - Z(j|(x)dx = Dp (17) 

x i 

where zj = desired negative altitude at the point x. 

Two or more quality criteria— Cases involving two or more quality criterion simultaneously 
(e.g., deviation minimum and fuel minimum) were also considered during the determination of 
optimal landing trajectories. The total performance index is then the sum of the squared and 
differently weighted single performance indices. 


RESULTS AND DISCUSSION 


The calculations of optimal landing flight paths for the augmentor wing STOL aircraft fall into 
three general categories. The first contains optimal input controls for straight line flight paths. The 
second category contains results obtained for optimum landing profiles on the basis of different 
individual criteria. These profiles are chiefly of theoretical interest. The third category contains 
more practical results. The control functions are more representative of those a pilot would choose 
to employ and the constraints are intended to represent those required in practice. This division of 
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results into theoretical and practical turns out to be somewhat arbitrary, because although the 
practical results begin with a trajectory with a complete set of constraints, the effects of removing 
some of them are then evaluated. 


Optimum Input Controls for Straight Line Flight Paths 

It is possible to determine optimal input controls even for straight line trajectories. Indeed, 
when fixing both the flight path angle and the speed of the aircraft, there is still one degree of 
freedom among the three control variables. 



Figure 12.- Fuel and noise minimum horizontal flight 
(5p = 75°). 


Horizontal flight paths— The aircraft is 
considered to be in equilibrium flying a straight 
horizontal course with its flaps deflected to 75°. 
Each control function has only a single 
parameter to be optimized for given quality 
criteria: pitch angle 0 = c(l,l), thrust angle 
a = c(2,l), and throttle angle 0 e = c(3,l). Results 
are shown in figure 1 2 where the abscissa gives 
the forward speed in knots, and the ordinate 
gives the controls 0, o, and 0 e for minimum fuel 
and noise, the state (a), the fuel required (Fp), 
and the average noise (Np) above 70 dB 
produced at an altitude of 700 ft (214 m) above 
ground for l.On.mi. (1.85 km). The solid lines 
show the values determined with a maximum 
angle of attack a = 8° ; the broken lines were 
obtained without this restriction. 

For a fixed speed, the indicated values 
correspond to both a horizontal minimum fuel 
flight from one given point to another, and to a 
horizontal minimum noise flight because when 
the altitude and the speed are fixed, noise 
depends only on the magnitude of the thrust. As 
an example of typical control settings, in 
horizontal equilibrium flight at V = 80 knots, 
figure 12 gives 0 = 3.19°, o = 18.58°, and 
0 e = 18.34°. 

Figure 1 2 also shows that a speed of about 
90 knots would result in a minimum fuel flight, 
and a speed of about 80 knots would result in a 
minimum noise flight. 
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Descent along a constant 7.5° glide slope- For a straight line descent along a given flight path at 
a given total speed, one degree of freedom remains for adjusting the three control variables (flap 
angle = 75°). Figure 13 shows these technically possible combinations of the control positions 0, o, 
and 0 e for 7 = -7.5° and V = 60 knots, with the resulting values of the consumed fuel and effective 
perceived noise when flying from an initial point 
(xj;zj) = (—5414 ft (-16 50 m); -700 ft (-2 14 m)) 
to the aim point (xf,zf) = (—100 ft (—31 m); 0). 

Theoretically it may be feasible to fly with a 
throttle angle of 0 e =lO° and corresponding 
0 = 14.76° and o = 34.63°. The resulting fuel 
consumption would be 41.8 lb (18.2 kg) and the 
average noise 98.7 EPNdB. 

Although this adjustment of the controls is 
advantageous with regard to the fuel and noise 
criteria, the safety constraints of d Q > 16° and 
a < 8° are not met. The a requirement is the 
more stringent and to meet it, one has to adjust 
input controls as follows: 0 e = 19.82°, 

0 = 0.42°, and o = 83.75° (fig. 13). As shown, 
the fuel requirement has now nearly doubled to 
Fp = 81.7 lb (35.6 kg) and the noise has 
increased to Np = 1 14.3 EPNdB, demonstrating the unfavorable effect on fuel consumption and 
noise propagation of the limits on the minimum throttle angle and maximum angle of attack. 

For a flight path angle of 7 = —7.5°, the fuel and noise optimal controls can be determined 
from figure 14, for any total speed 50 knots < V < 90 knots. Figure 14(a) shows the values when 
no restrictions are imposed on 0 e and a. Figure 14(b) shows how the controls have to be adjusted 
when the constraints are considered. It is apparent that for low speed the constraint on angle of 
attack a = 8° results in an increase of fuel and noise, and for high speed the constraint of a 
minimum throttle angle 0 e = 16° causes similar increases. 

If the constraints are taken into account, the least fuel and noise figures are obtained at 
approximately 90 knots. In the case of V = 60 knots, it is apparent that when the constraints are 
included the fuel usage increases from 41.8 (18.2) to 81.7 lb (35.6 kg) and the average noise from 
98.7 to 1 14.3 EPNdB (the figure shows the value above 70 dB). 
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Figure 13.- Inputs and states as a function of the throttle 
angle 8 e for 7 = -7.5° and V = 60 knots. 


Theoretical Optimum Landing Profiles 

For the determination of the controls and landing trajectories optimized with respect to noise, 
fuel, or time, the modified sine series (eqs. 9(a-c)) with eight coefficients for each input variable 
was chosen. The noise, fuel, or time during the whole trajectory from (xpzj) = (-1 .5 n.mi. 
(2.8 km);— 700 ft (214 m)) to (xpzf) = (0;0) was integrated using equations (1 1), (12), or (16b), as 
appropriate. 
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Velocity, knots 

(a) Without safety constraints . (b) With safety constraints on a and 0 e . 

Figure 14.— Fuel and noise minimal flight with 7 = -7.5° (5p = 75°). 

So that the aircraft might land at the desired touchdown point with the desired speeds 
(Vf;zf) = (55 knots; 3 ft/sec (0.9 m/sec)), the applicable performance criterion was augmented by a 
penalty function involving the weighted quadratic errors of the landing requirements. At the initial 
point the aircraft is assumed to fly in equilibrium with a horizontal speed of 80 knots and with the 
values of the controls indicated in figure 12. All the design limitations in table 1 were taken into 
consideration, but not the safety and comfort constraints in table 2. 

Noise minimum landing controls and flight paths - Figure 15 shows the best result of several 
optimization runs; there are numerous local optima since the problem is highly nonlinear. 

Looking at altitude as a function of the distance to the runway, it appears that the aircraft 
climbs from the initial 700 ft (214 m) to 875 ft (267 m). The thrust magnitude is highest during the 
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Figure 15.— Minimum noise landing trajectory with design and touchdown constraints only. 


first part of the flight. Thus it is profitable to fly as long as possible at the highest possible altitude, 
since the noise decreases as the distance to the noise source increases. During the middle part of the 
trajectory the thrust magnitude and corresponding noise level are very low (0 e going below 3°), the 
maximum perceived noise decreasing to 72 PNdB. At the end of the flight, however, the thrust has 
to be raised again to meet the touchdown requirements. For this reason, and partly because the 
distance of the aircraft to the ground is diminishing, the noise increases during the last part of the 
trajectory. The noise level drops sharply at x = -1 550 ft (-472 m) , the beginning of the clear zone in 
figure 1 1 . 

Note that the angle of attack is relatively high during the whole flight (up to 17°); it is more 
advantageous to get the necessary lift from a higher angle of attack than from a higher magnitude 
thrust. 

The average maximum effective perceived noise along this minimum noise flight trajectory is 
N = 96.7 EPNdB, the fuel consumption is Fp = 8 1 .5 lb (35 .5 kg), and the time needed to cover this 
1.5 n.mi. (2.8 km) flight path is Tp = 76.0 sec. 

Minimum fuel landing - The best minimun fuel controls and corresponding trajectories found 
are shown in figure 16. The thrust magnitude is extremely low at the beginning of the flight and the 
aircraft loses altitude quickly. For the remaining part of the trajectory the aircraft has to be pushed 
with a medium thrust to keep altitude and meet the final touchdown requirements. The average 
noise increases to Np = 103.4 EPNdB, while the fuel consumed is only Fp = 60.0 lb (26 kg) and the 
flight time decreases to Tp = 69.7 sec. 
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Figure 16.— Minimum fuel landing trajectory with design and touchdown constraints only. 

Minimum time landing - Whereas the noise minimal flight path showed at first a climb and the 
fuel minimal trajectory a dive, the time minimal flight path is nearly straight from the initial to the 
touchdown point (fig. 17). The thrust magnitude is very high over nearly the complete path to 
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Figure 17.— Minimum time landing trajectory with design and touchdown constraints only. 
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allow the aircraft to fly with the maximum permitted speed of V max = 95 knots (for 5p = 75°) at 
first and then to decelerate to reach an admissible touchdown speed. This reduction is performed by 
a thrust reversal (i.e., the thrust angle changes from approximately 22° to nearly 116°, the 
maximum allowable thrust angle). Another characteristic element of this trajectory worth 
mentioning is that a negative pitch angle is necessary during the entire flight. The greater thrust 
magnitude causes an increase both in speed and in lift. To allow the aircraft to lose altitude, one has 
to diminish the angle of attack by a negative pitch angle, particularly at landing when a throttle 
setting of 0 e = 34° combined with a thrust angle of nearly 1 16° causes high lift. 


For this minimum time flight the resulting 
average noise is Np = 108.9 EPNdB, the fuel 
required is Fp = 78.5 lb (34.2 kg) and the flight 
time is only Tp = 60.6 sec. In achieving these 
values, however, the optimization procedure for 
minimum time drove the system to a touchdown 
speed of 63.6 knots compared to the required 
speed of approximately 55 knots. This higher 
touchdown speed is nevertheless low enough to 
allow a proper landing on a 2000-ft (610 m) 
STOL runway (ref. 4). 

Recapitulation — Figure 18 shows the h/x 
trajectories obtained for noise, fuel, and time 
minimum flights, with corresponding values of 
the average noise produced, fuel required, and 
the time needed for these flights. 
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Figure 18.- 
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Landing trajectories without safety and 
comfort constraints. 


Optimum Landing Profiles With Realistic Constraints 

We have treated optimum flight paths considering only the design constraints of table 1 and 
touchdown conditions. In this section, profiles discussed include the constraints and conditions noted, 


as well as the safety and comfort constraints of 
from touchdown. Five landing paths are treated, 
a standard of comparison for the others. In the 
following trajectories, the effect on noise, fuel, 
and time of relaxing various constraints will be 
examined. Before examining these results, 
however, the height constraint and the control 
function structure must be explained. 

The height constraint will be explained in 
connection with the following discussion of the 
nominal landing flight. The nominal conditions 
are illustrated in figure 19. The aircraft starts in 
equilibrium at (x j ;zj) = (— 1 .5 n.mi. 
(-2.8 km)/-700 ft ( — 2 1 4 m)) with a speed of 


table 2 and a height limit at a point 350 ft (107 m) 
The first, satisfying all the constraints, will serve as 

I 
I 



Figure 19.- Desired landing flight trajectory. 
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Vj = 80 knots, 7 i = 0° and flap angle constant at 6p = 75° during the flight. It flies horizontally 
with a speed of 80 knots until its flight path crosses a fixed -7.5° glide slope. Then the total speed 
should be reduced to 60 knots, and the flight path has to follow as closely as possible the y = -7.5° 
glide slope until a point located at (x^jzjj) = (-350 ft (-107 m)/-35 ft (-1 1 m)) is reached. At this 
point the total airspeed should still be = 60 knots and the vertical speed is expected to be 
z^ = 13.22 ft/sec (4.03 m/sec), corresponding to the flight path angle of y = -7.5°, Finally, a flare is 
necessary, allowing a touchdown at (xfjzf) = (0;0) with desired speeds of Vf= 55 knots and 
Zf = 3 ft/sec (0.9 m/sec). To this nominal path we add the limitations on the control variables 
( a max> a min> l^lmax* ^e max > the constraints on the states (a max , V max ) and the 

constraints due to the comfort requirements of the passengers (|V| max , |0| max ) gi ven in tables 1 
and 2. 

A control-function structure was selected involving hyperbolic tangents as in equations (12a, b, 
and c) and illustrated in figure 10. These formulas sometimes involved four trajectory portions 
(three flight-path changes) and sometimes only three (two changes). The lengths and the centers of 
the switching intervals were taken to be the same for all three control functions. If a structure with 
three trajectory portions was chosen, two coefficients defined the center of the switching intervals, 
two more defined the length of these common switching intervals and 3X3 coefficients defined the 
amplitude of the steady state portions of the three input functions. Thus, the state optimization 
algorithm had to determine only the optimal values of 13 coefficients. When four trajectory 
portions were chosen, 3 + 3 + (3X4) = 18 coefficients had to be optimized. 

Fully constrained minimum noise landing flight path— When one sets all the limits in the 
optimization program and selects a quality criterion comprising a “noise minimum” plus a 
“deviation minimum” form, so that the flight remains as close as possible to the path specified on 
figure 19, the optimization procedure delivers the time histories of the inputs, states and outputs 
indicated in figure 20. (With the trajectory and speeds fixed, the minimum noise and minimum fuel 



- 1.40 - 1.20 - 1.00 -.80 -.60 -.40 -.20 0 

Distance to touchdown, nmi 


Figure 20.— Minimum noise landing controls and states involving all constraints. 
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controls are identical.) These results were obtained with input functions involving only three 
trajectory portions (horizontal 80 knot flight, -7.5° descent at 60 knots, and flare). Note that 
during the first part of the flight the aircraft is effectively flying with the desired horizontal speed 
(V = 80 knots, 7 = 0°), with application of the optimal controls (0 = 3.19, a = 18.58, 0 e = 18.34) 
which are the same as those shown in figure 12 . 

After a maneuver — performed so as to avoid violating the imposed technical and comfort 
constraints — the aircraft captures the second part of the trajectory, a flight path with y - -7.5° and 
V = 60 knots. Again the controls encountered (0 - 0.45, a = 83.74, 0 e = 19.81) correspond to the 
controls discussed in connection with optimal straight-line flight paths (see fig. 14(b)). Note that 
the angle of attack is 7.97°, very close to the limit. At the “flare initiation point” the aircraft is 
0.1 ft (0.03 m) lower than the required altitude of 35 ft (1 1 m) and the effective speeds at this 
point are V = 59.98 knots and z = 12.89 ft/sec (3.93 m/sec), which are nearly the desired speeds. 

The third and last part of this trajectory is the flare. With a nose up position, the aircraft 
touches down at Xf;zf = 0.04 ft (0.01 m)/~0.08 ft (-0.02 m) with the speeds of Vf = 57.18 knots 
and Zf = 2.98 ft/sec (0.91 m/sec), values that are again very close to the desired numbers. 

The average noise produced for this flight involving all constraints is Np = 1 10.0 EPNdB, the 
fuel consumed is Fp = 128 lb (56 kg) and the required time is Tp = 82.3 sec. 

Landing path with higher descent speed- A descent speed of 60 knots with a flight path angle 
of 7 = -7.5°, giving a vertical speed of z= 13.22 ft/sec (4.03 m/sec), is considered standard for a 
STOL aircraft landing trajectory. In figure 21, however, the optimal controls and corresponding 
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Figure 21.— Controls and states corresponding to a higher descent speed. 
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states are shown for a flight path angle of 7 = -7.5° with the higher descent speed of V = 80 knots 
and vertical speed of z = 17.62 ft/sec (5.37 m/sec). Here, the input control functions are divided in 
four portions: horizontal 80-knot flight, -7.5° flight path at 80 knots, switching to 60 knots, and 
flare. 


The resulting main differences compared to the basic fully constrained trajectory are: 


1 . Descent speed increased to about 80 knots 

2. Angle of attack during the descent much lower with a around 2° 

3. A momentary leveling off when the controls change to reduce the speed before touchdown 


The average noise resulting from the increased descent speed is slightly lower than for the basic 
trajectory with Np = 109.4 EPNdB. The fuel consumed decreases to Fp = 101.8 lb (44.4 kg) and 
the required time is smaller with Tp = 7 1 .0 sec. 

Landing path with a higher angle of attack— Vox adequate stall margin against vertical gusts it is 
reasonable not to permit descent with an angle of attack greater than 8 ° . This safety requirement 
has a strong influence on the noise and fuel characteristics, however (cf. figs. 20 and 22). The results 
shown in figure 22 required input functions involving only three trajectory portions. 
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Figure 22.- Controls and states corresponding to a landing with a higher maximum angle of attack. 

The controls and states in figure 22 are substantially the same as those in the basic fully 
constrained trajectory except for the angle of attack (and the associated controls), which increases 
to nearly 12°. The resulting noise, however, is 3.6 dB lower (Np = 106.4 EPNdB) the fuel is 
15.3 percent lower (Fp = 108.5 lb) (47.3 kg) and the time is practically the same (Tp = 82.5 sec). 
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Landing path with a higher descent speed, steeper flight path angle, and a lower throttle angle- 
in the example of figure 23, three constraints are relaxed: The descent speed is allowed to go up 
around 80 knots; the flight path angle, with y = —10°, is steeper than the recommended -7.5°; and 
the throttle angle is allowed to go nearly as low as 8°. 
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Figure 23 — Controls and states corresponding to a higher descent speed, a steeper flight path angle, and a lower 

throttle angle. 

For this trajectory, the input functions required four portions. The results shown in figure 23 
exhibit two main differences from the basic flight path (fig. 20): the aircraft momentarily levels off 
when changing controls to reduce speed before touchdown, and altitude at inception of flare is only 
33.0 ft (10 m),_rather than 35 ft (11 m). With the three constraints relaxed, Np is a very low 
102.8 EPNdB, Fp = 86.0 lb (37.5 kg) and Tp = 70.8 sec. 

Landing path with limited constraints — The final minimum noise landing path considered in 
this investigation was subject only to the constraints of table 1 and an initial horizontal speed of 
80 knots, which is maintained up to the point of intercepting the glide slope. The input functions 
were again divided into four portions. The results (fig. 24) differ significantly in several respects 
from those for the basic trajectory (fig. 20). Forward speed during descent is 80 knots rather than 
the nominal 60 knots, the flight path angle steepens to -11° instead of -7.5° and the angle of 
attack to 20° instead of 8°, and the throttle angle is reduced from 16° to 4.5°. Conditions at the 
nominal 35-ft (1 1 m) point are violated: Altitude is 73.4 ft (22.4 m) rather than 35 ft (1 1 m); the 
speed is 72.1 knots; and the sink rate is 22.1 ft/sec (6.74 m/sec), nearly twice the nominal value. 
Finally, the touchdown speed is also high — 70.2 knots instead of the nominal 55 knots — which 
still is satisfactory for landing on a 2000 ft (610 m) runway (ref. 4). 

By integrating the noise developed along this bold landing trajectory and taking its average, 
one obtains Np = 95.6 EPNdB, 14.4 dB lower than the basic trajectory value of 
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Np = llO.OEPNdB; 1 the fuel consumed Fp becomes 70.6 lb (30.8 kg) (44.8 percent less than the 
basic trajectory value) and the flight time Tp is 69.1 sec, minus 13.2 sec lower than the basic value. 
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Figure 24.- Controls and states when considering limited constraints. 

Recapitulation - Figure 25 summarizes the results of the five cases described in this section. The 
main characteristics of the different landing flights are specified on the left and the resulting average 
noise produced, the fuel consumed, the time required, and the change in noise, fuel, and time from 
the basic trajectory are shown on the right. 


Computational Requirements for Optimization of the Flight Paths 

The entire Fortran IV computer program consists of four decks. The first, or main, deck, 
which comprises 80 cards, permits card reading, definition of the input data, and data printout. The 
data may be: constants, the variable characteristics of the aircraft, given initial and desired final 
values of the trajectory to be optimized, initial guesses for the control function coefficients, and 
others. 

The computational flow is then directed to the second program, the static optimization 
algorithm EXTREM, described in the appendix, with 100 Fortran statements. 


'This average noise produced is lower than the value of Np= 96.7 EPNdB encountered during the determination of noise 
minimal landing flights with the use of sine-series input functions (cf. fig. 15). In that investigation the safety and comfort constraints 
were not considered, but the touchdown requirements nevertheless were met: namely, the speed was held at a low level with 
Vf = 54.5 knots, whereas, here the touchdown speed is allowed to be very much higher with Vf = 70.2 knots. 
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After each specification of a new set of 
input coefficients, the subroutine EXTREM calls 
a third program, which evaluates the appropriate 
performance index. To do this the trajectory has 
to be integrated from the initial to the final 
point. The noise, fuel, time, and flight path 
deviations are evaluated at each point along the 
trajectory and summed. The actual performance 
index is one or a weighted combination of these 
quality criteria. This Fortran subroutine permits 
the selection of the Adams-Moulton, 
Runge-Kutta, or Euler methods for integration 
and allows the evolution of the inputs, states, 
and outputs to be printed out and stored on 
disk. It comprises nearly 200 cards. 

The final program is called by the third 
program at each integration step. Here the input 
values of 8, o, and 0 e are determined from 
equations (9a-c) or (lOa-c). It also evaluates the 
lift and drag coefficients (figs. 9(a,b)) 
corresponding to the current state of the aircraft 
and delivers the derivatives u, w, x, and z 
(eqs. (4a-d)) needed for the flight path 
integration. This fourth program, which includes 
1 1 different structures for input functions, 
consists of 120 Fortran statements. 
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Figure 25.- Effects of constraints on noise, fuel, and time 
characteristics of_a landing trajectory. (In the last four 
cases, values of Np, Fp, and Tp given are differences 
from those of the basic trajectory.) 


Using single precision, the entire program of the four decks, with a total of about 500 cards, 
requires a storage capacity of 4220 hexadecimal bytes on the IBM 360/67 computer Operating 
System (OS). The storage requirements associated with the number of coefficients to be optimized 
is K(K+5) words where K is the number of coefficients. 


The computation time necessary for one optimization run depends on the complexity of the 
problem, the number and type of constraints involved, the number of coefficients to be optimized, 
and the choice of initial values of these coefficients. Disregarding its dependence on the structure of 
the control functions and on the choice of initial values, the time of computation increases roughly 
as the square of the number of input coefficients to be optimized. 


The problem of the optimum fuel flight path with design constraints only (fig. 16) can be used 
as an example of the time required for finding an optimum. The following initial values were used 
for coefficients of the functions given in equation (9): 


c( 1 , 1 ) = 16 ; c( 1 ,2) = 2.0 ; 

c(2,l) = 26 ; c(2,2) = 3.0 ; 

c(3,l) = 20; c(3,2) = 18; 


c( 1 ,3) = .* • • c(l,8) = 0 


c(2,3) = • • • c(2,8) = 0 


c(3,3) = • • • c(3,8) = 0 
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The initial search step size chosen for all 24 coefficients was |DC| = 1 . Convergence with reasonable 
precision to the final trajectories shown in figure 16 took 7 stages (see appendix) or about 600 
trials. One trial corresponds to the integration of one complete trajectory, or one evaluation of the 
performance index, which is equivalent. This optimization run requires about 5 min using Euler 
integration with a step size of At = 1 .0 sec. 


CONCLUSIONS 


An engineering dynamic optimization method was described that requires access only to the 
inputs and the output of the system whose performance is to be optimized. The method consists in 
transforming the dynamic problem into a static one by representing the input controls as analytical 
functions whose coefficients have to be determined by a parameter optimization procedure that 
minimizes or maximizes certain quality criteria adjoined to the system. 

The method was applied to the determination of open loop commands for a STOL aircraft in 
the landing phase. The criteria of noise, fuel, time, and flight path deviation were considered as 
performance indexes. 

The problem required the simultaneous time history optimization of three control variables — 
pitch angle, thrust angle, and throttle angle — under the influence of input, state, and terminal 
constraints. Three types of results were obtained: 

1 . Optimal controls for straight-line flight paths 

2. Optimal controls and corresponding landing flight paths where only the design constraints 
of the aircraft and terminal conditions were considered 

3. Optimal landing controls and trajectories of practical interest, in the sense that safety, 
comfort, and flight path constraints were added to the design constraints. 

The first type of results indicated that, for a given flight path angle and speed, one should 
adjust the throttle angle to the lowest admissible value satisfying both the constraints of minimum 
throttle angle and maximum angle of attack, for both minimum noise and minimum fuel criteria. 
The trajectories obtained with only design and touchdown constraints show that for a noise 
minimum landing flight it is reasonable to fly at first as high as possible for as long as possible, then 
to descend with a minimum thrust amplitude and a maximum angle of attack, and finally to 
increase thrust for touchdown. 

If fuel consumption is the optimization criterion, then the flight path trajectory shows an 
initial rapid descent requiring very little fuel. For the middle and last part of the flight a medium 
fuel rate is necessary to hold the altitude and meet the final conditions. 

In contrast, the time minimum flight path trajectory displays a nearly straight-line shape. 
During the complete flight the thrust magnitude is high at first to accelerate the aircraft to the 
allowed speed limit and then to decelerate it, with thrust reversal, to an admissible touchdown 
speed. Another characteristic item of this flight is a negative pitch angle for the whole trajectory. 
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For trajectories involving safety, comfort, and flight path constraints, in addition to the 
touchdown and design constraints, it was concluded that: to decrease the average effective noise 
during a landing, the descent speed, angle of attack, glide slope, and touchdown speed should be 
held at the highest possible or highest admissible values. Because of the constraints imposed, the 
fuel and time required also tend to be reduced. 

Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., 94035, June 1, 1972 


27 



APPENDIX A 


FINDING AN EXTREMUM OF A BOUNDED MULTIVARIABLE FUNCTION 
WITHOUT DETERMINATION OF THE DERIVATIVES 


The static optimization algorithm is explained here in more detail. (The corresponding 
Fortran-Subroutine EXTREM is displayed at the end of this appendix.) This algorithm is best suited 
for the search of a local extremum (maximum or minimum) of a multivariable function, the gradient 
of which is impossible or difficult to obtain directly (e.g., functional optimization or parameter 
optimization of a complex system). The function, which may be analytical or computed indirectly, 
can include any number of independent variables and any number of inequality constraints on them 
or functions of them. It is possible to change the boundaries during the optimization procedure. 


METHOD 


Basically, the algorithm for the search of the extremum of a bounded multivariable function 
has to perform the following tasks: (1) choose the directions of search, (2) determine the optimum 
along a line, (3) define the size of each search step, and (4) consider the boundaries limiting the 
search procedure. 


Choice of Search Direction 

The first main line (optimal search direction) is given by the user in form of an array DX(K), K 
being the number of variables. This array contains K properly scaled increments defining the initial 
step sizes along each variable about the guessed initial point Xj, which is defined as a vector in an 
array X(K). Along this first line the approximately extremal point Xj +1 is determined. After this 


Secondary direction 
of first stage 



Figure 26.— Determination of the main search direction in 
the case of a two variable problem. 


first iteration a second line, going through Xj+! 
and being orthogonal to the main direction, is 
defined by a Gram-Schmidt orthogonalization 
process (fig. 26). Again the near-extremal point 
Xj +2 along this secondary direction is 
determined (second iteration) and a third line is 
defined, if there are more than two coefficients, 
that passes through Xj +2 and is orthogonal to 
the two first directions. 

This procedure is repeated until the 
extremal point Xj+k along the Kth line, which is 
orthogonal to all previous directions, has been 
determined. Thus the first stage (K iterations) is 
accomplished. The new main direction for the 
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following stage is now given by the line going through the initial point Xj and the extremal point 
Xj+k of the last iteration of the previous stage. The procedure just described for the determination 
of the main and secondary search directions may begin again by setting Xj = Xj+fc. 


Determination of the Optimum Along a Line 


To start, the function values, Fj at X, - DX, F 2 at X^ and F 3 at 5^ + are evaluated 
(search steps) and a parabolic extrapolation (or interpolation) is performed to find the extremum of 
a fictitious parabola going through the three defined points (extrapolation or interpolation step) 


^i+i Xj + 


DX 


IF, 


2F 2 + F 3 


f 3 -F, 

2M m 


where 


Mm = + 1 for the search of a maximum 


= -1 for the search of a minimum 


This equation is obtained by solving for the 
coefficients of a parabola at the three given 
values, then determining the position of Xj+x at 
the extremum of the parabola. In this formula 
the absolute value of (F, — 2F 2 + F 3 ) is taken 
so that the new point Xj +1 will be closer to the 
expected extremum than the old point Xj even 
if inflection points are located between the 
identified part of the curve and the extremum 
(see fig. 27). 



Figure 27.— One-dimensional case when an inflection point 
is located between the identified part of the curve and 
the extremum sought. 


The algorithm limits the progression of the 
new point Xj+i with respect to a maximum of 
20 step sizes M (e.g., if IF, -2F 2 + F 3 | = 0). 

The function value F 0 pj at the new near-optimal point Xj+j along the line is evaluated, and these 
coordinates are taken as bases for the next iteration, even if F 0 pj is worse than F 2 . However, if an 
undesirable difference between F 0 pt and F 2 is greater than four times the absolute value of the 
difference between the function values at the end of the last two iterations, the progression of the 
new point Xj+, is divided by two and another function value is evaluated. The same procedure is 
repeated until the above criterion for the transition to the next iteration along another line holds. 
Note that in some cases where the new point is allowed to have a slightly worse function value than 
the old point, the operating point leaves the bottoms of narrow valleys or the tops of sharp ridges, 
which is advantageous when these valleys or ridges present sharp comers. 


Definition of Step Size 

The different step sizes along each direction for the first stage are given by the array DX(K). If 
during any iteration the distance between the new point and the old is smaller than one quarter of 
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the current step size along this line, the current step size is divided by 4. If the new point is more 
than 20 times the current step size away from the old point, then the step size will be multiplied by 2. 

By this procedure the search steps are allowed to decrease near the optimum or within tight 
curves and to increase again once these curves are passed. 


Consideration of Boundaries 

Before any function value is computed, the program given by the user is checked to determine 
if the current point lies within the permitted area. If this is not the case, the control is returned to 
the program EXTREM, which will deliver another point that satisfies the constraints. To compute a 



new point, if the previous one was outside the 
permissible area, the algorithm discerns two 
cases: (1) if the previous trial has been a search 
step as outlined in the preceding section, the 
step size is divided by 4 and the new point Xj+j 
will be placed at the other side of Xj away from 
the boundary (fig. 28, line A) and (2) if the 
previous trial has been an extrapolation or 
interpolation step, the progression of the new 
point Xj+! in the direction of the boundary will 
be divided by 10 (fig. 28, line B). 


Restrictions 

The guessed initial vector, given by the user, should not lie outside the defined boundaries. 
Moreover, even though one may shift the limits defining the valid region while the operating point is 
converging to the extremum, current values of the variables (arguments) must not lie beyond the 
boundaries. 

The stopping conditions (on the function value, the arguments, and the number of stages 
allowed to reach the optimum) should be sufficiently severe (small values for DFMAX and 
DXMAX, and a large number for LMAX) to avoid false convergence (e.g., stopping at sharp edges). 
This severeness is reasonable, especially when sharp edges exist. One must allow the step size to 
diminish sufficiently so that the optimization procedure, instead of stopping at the barriers, may 
follow the contour of the edges. 


APPLICATION OF THE ALGORITHM 


Calling Sequence 


EXTERNAL 

FN 

DIMENSI0N X(...), DX(...), S( ) 

CALL EXTREM (FN, K, X, DX, S, DFMAX, DXMAX, LMAX, F0PT, IW) 
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with FN - name of subroutine supplied by user to determine if 
current arguments lie within the specified boundaries and, 
if yes, to evaluate the corresponding function value. 

K - positive number whose magnitude is the number of 
independent variables of the function. 

X - one-dimensional array whose K elements contain the initial 
guess arguments of the K independent variables. At the end 
of the optimization they deliver the final values of those 
arguments. 

DX - one-dimensional array whose K elements define the K 
initial step sizes along each variable about the initially 
guessed point. 

S - two-dimensional array defining the supplementary working 
space needed; it should be dimensioned (K, K + 3). 

DFMAX - stopping condition on function variation; the optimization 
procedure stops if during the last stage the variation of the 
function value was smaller than DFMAX. 

DXMAX - stopping condition on arguments; the optimization 
procedure stops if during the last stage the absolute 
variation of the argument vector was smaller than 
DXMAX. 

LMAX - stopping condition on number of stages; the optimization 
procedure stops if the current number of stages is equal 
|LMAX|. The sign of LMAX indicates if a maximum is 
sought (positive sign) or a minimum (negative sign). 


F0PT - function value at the end of the optimization procedure. 

IW - Writing instructions: 

±1 all outputs suppressed (e.g., results transferred to the 
user’s main program). 

±2 final outputs only 

±3 outputs at the end of each stage 

The sign of this instruction indicates if boundaries are involved 
(positive sign) or not (negative sign). Note that the optimization 
procedure stops if at least one of the three stopping conditions holds. 



Form of the Output 


STAGE NUMBER ... TRIAL NUMBER .... DL = 

FU = AR(1)= DS(1) = 

AR(2) = DS(2) = 


AR(K) = DS(K) = 


where 

DL 

FU 

AR(1) . . . AR(K) 
DS(1) . . . DS(K) 


= magnitude of the argument vector variation during the last stage 
= current function value 
= current argument values (variables) 

= step sizes in the main direction S( 1 ) and in the orthogonal secondary 
directions S(2) . . . S(K) 


The user must supply a subroutine for the determination of the function values. Any 
inequality constraints on the arguments or on functions of the arguments may be introduced in this 
subroutine by setting LI = LI + 1 (defined below) in the case that one or several of the current 
arguments lie on the wrong side of the specified boundaries. In this case a special instruction sends 
the flow back to the program EXTREM and the function value is not computed. New arguments are 
then determined by the program EXTREM, which do not violate the boundaries. This subroutine 
FN may have the following form: 

SUBROUTINE FN(AR, F, LI, N) 

DIMENSI0N AR(1) 

IF (LOGICAL EXPRESSION) LI - LI + 1 I necessary on j y j n case that boundaries are involved 
IF (LI.GT.l) RETURN I 

F = ... 

N = N+l 

RETURN 

END 

with: 

LOGICAL EXPRESSION - should be true if one or several of the arguments lie in the forbidden 

area 

F - function value corresponding to the current arguments 
N - current number of trials (function evaluations) 

The constraints may be changed at each trial; as noted, they can be a function of the current 
arguments AR(K) or even of the current function value F. 
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Storage Requirements 


The program contains built-in variable dimensions so that its working space is defined by the 
calling program; the total vector argument storage needed (X, DX, and S vectors) is K*(K+5) for K 
independent variables and for any number of inequality constraints. The program itself comprises 
fewer than 100 instructions and takes 15A4 hexadecimal bytes of storage (double precision) for the 
Fortran G compiler of the OS. 


EXAMPLES 


The extrema of both analytic and computed functions have been determined with the 
described method. The optima-search results of five analytic functions are shown (computation in 
double precision) here. The first is Rosenbrock’s curved valley (ref. 1 1 ): 

F = 1 00(Xj 2 — X 2 ) 2 + ( 1 — X, ) 2 ( A1 ) 


where 

True min: Xj = X 2 = 1 , F = 0 
Initial values: X l =-1.2, X 2 = 1 
Initial step sizes: DX! = 1, DX 2 = 1 

State after 102 trials (function evaluations): X t = 0.9858, X 2 = 0.9718, F = 2.006X10 -4 . 

The second is a test function suggested by Aubrun (ref. 5), a steep, narrow, winding ridge with 
a horizontal inflection point along the top of this ridge: 

F = -100(Xj 3 + X, 2 -X, — X 2 ) 2 -(X, 3 -l) 2 < A2 > 


where 

True min: X, = X 2 = 1 ; horizontal inflection point Xj = X 2 = 0 
Initial point: X t = -2.2, X 2 = -2 

Initial step sizes: DX! = DX 2 = 1 

State after 368 trials: X, = 1.0011, X 2 = 1.0046, F= 1.2237X10 -5 . 

The third is a linear function with discontinuities which cannot be approximated by a 
quadratic — given by R. F. Wheeling (ref. 12): 

F = — 3 |X, | — |X 2 1 (A3) 

where 

Max: X, = X 2 = 0 

Initial point: Xj = X 2 = 10 

Initial step sizes: DXi = DX 2 = 1 

State after 1 29 trials: Xi = -7.8899X 1 0 -5 , X 2 = -2.5 1 47X 1 0" 4 , F = -4.88 1 7X 1 0 -4 



The fourth is Rosenbrock’s parcel problem where the maximum lies on three boundaries 
(ref. 11 ): 


F = X t X 2 X 3 with boundaries 


0 <Xj <20 

0 < X 2 < 1 1 

0 < Xj + 2X 2 + 2X 3 < 72 , 


(A4) 


where 

Correct result: X! = 20; X 2 = 1 1 ; X 3 = 15, F = 3,300 
Initial point: X 1 = X 2 = X 3 = 10 
Initial step sizes: DX, = DX 2 = DX 3 = 1 

State after 484 trials: X, = 19.9999, X 2 = 10.9997, X 3 = 15.0001, F = 3,299.953. 

The fifth is Powell’s fourth-power function, which cannot be approximated by a quadratic 
near the minimum (ref. 13): 

F = (X t + 10X 2 ) 2 + 5(X 3 -X 4) 2 + (X 2 — 2X 3 ) 4 + 10(Xj -X ^,) 4 (A5) 

where 

Min X, = X 2 = X 3 = X 4 = 0 

Initial point: Xj = 3;X 2 = -1 ; X 3 = 0; X 4 = 1 

Initial step sizes: DXj = DX 2 = DX 3 = DX 4 = 1 

State after 181 trials: X t = -9.7493X10 -3 , X 2 = 9.7492X1 0~\ X 3 = -3.7262X 10 -3 , 
X 4 =-3.7280X1 (T 3 ,F = 1 .8205X10”®. 

Further reduction of F is very slow. 

The program EXTREM has also been used to search the optimum of some bounded 
functionals (computed functions), such as the optimal trajectories of rockets and the STOL aircraft 
described in this report. The procedure was tried successfully on problems involving up to 30 
parameters to be optimized. 


FORTRAN PROGRAM 


Despite its simplicity and relatively small size, the Fortran program given below finds the local 
minimum or maximum nearest a given initial point of any constrained analytical or computed 
multivariable function or functional. This algorithm is believed to be more efficient (in terms of the 
number of times the function must be evaluated to reach the extremum with a given precision) than 
other procedures described in the literature (ref. 1 4). 
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SUBROUTINE EXTREM < F , K, X, DX , S , DFMAX , DXMAX,LMAX, FOPT , IW ) 

HGJ 

01 

c 

FINDING AN EXTREMUM OF A BOUNDED MULTIVARIABLE FUNCTION 

HGJ 

02 

c 

WITHOUT DETERMINATION OF THE DERIVATIVES 

HGJ 

03 

c 

(SINGLE PRECISION PROGRAM) 

HGJ 

04 



DIMENSION X(1),DX(1),S(K,1> 

HGJ 

05 



L=0 

HGJ 

06 



LI = 1 

HGJ 

07 



N=0 

HGJ 

08 



DO 1 1=1, K 

HGJ 

09 



S ( I » 1 ) =X ( I ) 

HGJ 

10 


1 

S< 1,2 ) = X ( I ) -DX ( I ) 

HGJ 

11 



CALL F(X,F2,LI,N) 

HGJ 

12 



FE=F2 

HGJ 

13 



IF(LI.GT.1)WRITE(6,2) 

HGJ 

14 


2 

FORMAT( IX, • INITIAL ARGUMENTS OUTSIDE BOUNDARIES') 

HGJ 

15 


3 

IF(KC.GE.K.0R.KC.LT.0.0R.L.EQ.0)KC=0 

HGJ 

16 



KC=KC+ 1 

HGJ 

17 



S< 1,3)=0. 

HGJ 

18 



DO 4 1=1, K 

HGJ 

19 



S ( I ,4) =S ( I » 1 )— S ( 1,2) 

HGJ 

20 


4 

S(1,3)=S(1,3)+S(I,4)**2 

HGJ 

21 



S ( 1 , 3 ) = SQRT(S(1,3>) 

HGJ 

S22 

c 


S( 1,3 )=DSQRT(S( 1,3) ) 

HGJ 

D22 



I F ( IABS( IW) .GE.3)WRITE( 6,23)L,N,S(1,3) »F2»( I ,X( I) ,1 ,DX( I ) , 1 = 1 ,K ) 

HGJ 

23 



IF(L.GE.IABS(LMAX).OR.S( 1, 3) .LT.DXMAX .OR. ABS< FF-FOPT) 

HGJ 

S24 

c 


IF( L.GE. IABS(LMAX) .OR.S< 1,3) .LT .DXMAX .OR .DABS( FF-FDPT) 

HGJ 

D24 



l.LT.DFMAX.AND.L.GT.O.OR. (LI.GT.l.AND.L.EO.O) ) GOT 022 

HGJ 

25 



I F ( K .E Q. 1 ) G0T09 

HGJ 

26 



DO 8 J=2 , K 

HGJ 

27 



KD=— 2+J+KC 

HGJ 

28 



IFIKD.GT.K )KD=KD-K 

HGJ 

29 



S( J,3)=0. 

HGJ 

30 



DO 7 I = 1 , K 

HGJ 

31 



S ( I » J+3) =0. 

HGJ 

32 



I F ( I.EO.KD)S( I, J+3)=S( 1,3) 

HGJ 

33 



JM= J— 1 

HGJ 

34 



DO 6 J K= 1 , JM 

HGJ 

35 


6 

S( I ,J+3)=S( I ,J+3)-S(KD,JK+3)*S( 1,3)/S( JK,3)*S( I , JK+3 ) / S ( JK ,3 ) 

HGJ 

36 


7 

S( J,3 )=S( J,3)+S(I,J+3)**2 

HGJ 

37 



S ( J ,3 ) = SORT ( S ( J ,3 ) ) 

HGJ 

S38 

c 


S( J,3 ) =OSQRT ( S ( J , 3 ) ) 

HGJ 

D38 



I F ( S ( J,3).LT.1.D-30)G0T03 

HGJ 

39 


8 

CONTINUE 

HGJ 

40 


9 

DO 10 1=1, K 

HGJ 

41 


10 

S( 1,2 )=S( 1,1) 

HGJ 

42 



L = L+ 1 

HGJ 

43 


35 ' 



FF=FOPT 
DO 21 M=1,K 
DO 11 I = 1 » K 

11 S( I »M+3 ) = S ( I ,M+3)/S(M,3 )*DX ( M ) 

IF ( I W .GT .0 ) L I =3 

12 IF ( IW.GT.0)LI=LI-1 
L J=L I 

DO 13 1= 1 »K 
X ( I ) =S ( I , 1 )~S ( I , M+3 ) 

13 S(I,M+3)=S(I,1 )— X( I ) 

CALL F<X,Fl,LItN> 

B0=1. 

14 DO 15 I = 1 1 K 

X( I)=S( Itl)+S( I » M+3) /BO 

15 S(ItM+3)=X(I)-S(I,l) 

I F ( ABS(BO) .GT. 1 .1 ) GOTO 20 
C IF(DABS(B0).GT.1.1) G0T020 

CALL F(X,F3,LJ,N) 

IFILI+LJ.EQ.4)G0T012 
IF(LJ.GT.2)B0=— 4. 

IF(LI.GT.2)B0=+4. 

IFILI.GT .2 .OR .L J .GT. 2) GOT 014 

16 ST=0. 

DO 18 1 = 1 ,K 
XIIIsSlUl) 

I F ( ABS( S( I »M+3) ) . LT . 1 . 6-30) GOTO 18 
C I F ( DA BS ( S ( I t M+3 ) ) . L T. 1 . D— 30 ) GOTO 18 

S ( I »M+ 3 ) =S ( I »M+3)/LI 

I F ( ABS ( 2 • *F2— F 1-F3 ) . LT. 1.6-30) GOTO 18 
C I F ( DABS ( 2 .*F2— Fl-F 3).LT .l.D-30)G0T018 

X( I >=S( I t 1 )+S < I ,M+3 )/ ABS(F1-2.*F2+F3)*(F3-F1)/ISIGN(2,LMAX) 
C X( I )=S( I »1)+S< I ,M+3 )/DABS(Fl-2.*F2+F3)<=(F3-Fl)/ISIGN(2,LMAX) 

18 ST=ST+(X< I )-S( 1,1) )**2 

IF( 16.*ST.LT.DX(M)**2)DX(M)=DX<M)/4. 

IF( ST .LT.400. -DX ( M ) **2. AND. ABS( 2. *F2-F1-F3).GE. 1.6-30) G0T020 
C IF(ST.LT.400.*DX(M)* C *2.AND.DABS(2.*F2-F1-F3).GE.1.D-30)GOT020 

DO 19 1=1 ,K 

I F ( ABS ( S ( I ,M+3) ).LT. 1.6-30) GOTO 19 
C IF(DABS(S( I»M+3) ) .LT.1.D-30)G0T019 

X(I)=S(I,1)+ S IGN ( S( I ,M+3 ) ♦ ( F3— F 1 ) /S ( I , M+3 ) ) * I S IGN( 20 »LMAX ) 

C X(I)»S( I, 1 )+DSIGN(SI I,M+3), (F3-F1)/S( I, M+3) )* IS IGN ( 20 ,LMAX ) 

19 CONTINUE 

DX ( M ) =DX ( M ) *2 . 

20 L I =+ 1 
BO=-BO 

I F { ABS(BO) .GT.1.1 )0X(M)=DX(M)/3. 
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C IF(DABS(80).GT. l.l)DX(M)=DX(M)/3. 

CALL F(X,FOPT,LI ,N) 

I F ( L I .GT .1 ) L I = 10 

I F ( ISIGN(1»LMAX)*( F0PT-F2 ) .LT ABS( FE-F2 ) *4 . . AND .L I .NE . 10 )L I =2 
C IF(ISIGN(1,LMAX)*(F0PT-F2).LT.-DABS(FE-F2)*4.. AND.L I.NE.10)LI=2 

IF(LI.GT.l. AND. AB S< BO > .GT . 1 . 1 ) G0T014 
C I F ( L I .GT. 1 . AND. DABS ! BO ) .GT. 1 • 1 ) G0T014 

IF(LI.GT.1)G0T016 
FE = F2 
F2=F0PT 
DO 21 1=1, K 

21 S(I,1)=X(I) 

G0T03 

22 I F ( I ABS( IW) .EQ. 2) WRITE! 6,23)L,N, S( 1,3) , F2 , ( I , X (I) , I , DX ( I ) , 1= 1 ,K ) 

23 F0RMAT(//1X,9HSTAGE NO. , 13, 10X,9HTRIAL NO. , 16, 12X,3HDL=E15.8,/1X, 
13HFU=E 15.8,/(23X*3HAR(,I2»2H)=»E15.8»5X»3HDS(»I2,2H)=»E15.8) ) 

RETURN 

C THE LETTER D BEFORE THE IDENTIFICATION NUMBER MEANS DOUBLE PRECISION* 
C THE LETTER S SINGLE PRECISION. 

END 
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conducted so as to contribute ... to the expansion of human knowl- 
edge of pheno.mena in the atmosphere and space. The Administration 
shall provide for the ividest practicable and appropriate dissemination 
of information concerning its activities and the results thereof.” 

— National Aeronautics and Space Act of 1958 


NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 

TECHNICAL REPORTS: Scientific and 
technical information considered important, 
complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 

Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. 

CONTRACTOR REPORTS: Scientific and 
technical information generated under a NASA 
contract or grant and considered an important 
contribution to existing knowledge. 

Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION OFFICE 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include conference proceedings, 
monographs, data compilations, handbooks, 
sourcebooks, and special bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may be of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs, 
Technology Utilization Reports and 
Technology Surveys. 


